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This  report  is  a  chronological  documentation  of  the  research  progress  made  by  Martijn 
Kolloffel  throughout  the  Fall  Semester  2003.  The  research  focuses  on  the  use  of  cumulant 
functions  in  queueing  theory  and  stochastic  modeling.  This  report  is  a  result  of  a  DOD 
research  grant  proposed  by  Dr.  Timothy.  I.  Matis,  with  the  purpose  to  engage  undergraduate 
students  from  New  Mexico  State  University  in  research  in  stochastic  models.  My  progress 
was  monitored,  evaluated,  and  documented  through  the  participation  in  the  undergraduate 
research  course  IE  400.  The  investigation  into  the  application  of  cumulant  functions  in 
stochastic  modeling  is  a  continuation  of  research  activities  in  the  spring  semester  of  2003. 

The  usefulness  of  cumulant  based  analysis  methods  is  researched  by  the  formulation  of  a 
suitable  model.  The  first  goal  in  this  semesters  research  endeavor  is  the  definition  of  a  model 
that  is  relevant  to  military  applications.  After  the  model  is  completely  defined,  we  use  a 
cumulant  derivation  procedure  to  find  an  approximation  to  the  measure  of  interest.  We  then 
validate  our  solution  by  comparing  it  to  a  model  simulation.  After  validation  we  strive  to 
expand  the  scale  of  the  initial  model,  which  can  show  the  mathematical  tractability  of  this 
procedure  for  large-scale  systems.  We  also  intend  to  expand  the  model  by  using  phase  type 
distributions,  which  allows  us  to  model  most  probability  distributions. 

During  the  summer  of  2003  some  ideas  for  formulating  a  model  suitable  for  cumulant-based 
modeling  had  been  generated  in  relation  to  a  NASA  research  grant,  that  focused  on  the 
reliability  of  components  that  are  used  in  space  technology.  We  think  that  the  stochastic 
analysis  of  component  reliability  in  military  applications  can  also  be  relevant  and  useful  to 
the  Department  of  Defense.  Optimizing  maintenance  schedules  of  military  equipment  can 
prevent  critical  system  failure  and  improve  system  availability. 

The  first  papers  in  my  literature  research  include  the  “Analysis  of  equipment  availability 
under  varying  corrective  maintenance  models”  by  Cassady  and  Iyoob.  [1],  and  a  paper  by  M. 
Kijima:  “Some  results  for  repairable  systems  with  general  repair”  [2].  Reading  these  papers 
resulted  in  thinking  about  availability  or  reliability  of  systems.  The  initial  model  we 
formulate  is  a  gearbox  in  which  the  rate  of  failure  is  dependent  on  the  state  of  the  individual 
variables  in  the  gearbox.  These  variables  could  be  the  temperature,  pressure,  or  the  number 
of  particles  in  the  box.  We  assumed  that  these  variables  of  the  gearbox  are  in  a  state  of  no, 
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light,  moderate,  or  severe  wear.  We  are  not  completely  satisfied  with  our  formulation  and 
foresee  some  problems  when  applying  the  cumulative  procedure  to  this  model,  which 
inspires  me  to  step  up  the  library  research. 

I.J.  Rehmert  generates  my  interest  in  shock  models  after  reading  his  dissertation:  “Time 
dependent  availability  analysis  for  the  Quasi-Renewal  process”.  Each  shock  arriving  to  the 
system  causes  component  failure  that  is  a  function  of  the  total  accumulated  damage  from 
previous  shocks.  The  components  are  monitored  at  discrete  points  in  time,  which  can  be 
modeled  with  a  discrete  aging  process. 

This  leads  us  to  investigate  compartmentalized  wear  models  involving  crack  propagation.  We 
formulate  a  model  in  which  cracks  initiate  discretely  according  to  a  Poisson  distribution,  and 
propagate  according  to  a  continues  non-linear  function.  The  damage  at  time  t  can  be 
represented  as  an  integral  of  an  analytical  function,  the  “crack  growth  expression”.  I  tried 
viewing  crack  propagation  as  a  deteriorating  system,  but  once  again  encountered  an 
expansion  of  the  model  that  was  unwanted  at  this  point  in  time.  In  order  for  me  to  be  able  to 
construct  an  analytical  and  simulated  model  within  the  timeframe  given  I  opted  for  a  more 
simplistic  model,  which  would  stress  the  advantages  of  using  cumulative  derivation.  This 
brought  us  right  back  at  thinking  about  component  reliability. 

We  formulate  a  model  in  which  a  single  component  operates  in  an  environment  of  non-fatal 
repairable  shocks,  or  impacts.  These  shocks  represent  system  stresses  or  other  components 
failing  in  the  system.  The  shocks  will  cause  repairable  and  cumulative  damage  to  the 
component.  The  arrival  of  shocks  is  according  to  a  Poisson  distribution.  The  failure  of  the 
component  is  a  result  of  both  the  un-repaired  or  current  and  the  cumulative  shocks  at  time  t. 
We  can  write  an  intensity  function  corresponding  to  component  failure  that  is  a  nonlinear 
function  of  the  current  and  cumulative  number  of  shocks.  We  can  now  use  the  cumulant 
derivation  method  to  obtain  an  approximation  to  the  probability  that  the  component  is 
operational  at  time  t.  An  extended  explanation  of  the  model  described  above  can  be  found  in 
the  appendix. 

To  great  relief  we  finally  formulate  a  model  that  seems  to  fit  the  scope  of  my  research. 
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The  next  step  is  the  creation  of  numerical  results  based  on  our  analytical  model.  Mathematica 
software,  and  a  program  for  the  Cumulant  derivation  procedure,  which  Dr.  T.I.  Matis  has 
developed,  helps  us  solve  our  partial  differential  equations.  We  define  several  intensity 
functions  and  evaluate  the  availability  of  the  component  at  time  t,  letting  t  increase  discretely 
from  0  to  50.  The  results  are  imported  in  excel  to  make  a  graphical  comparison  with  the 
simulated  results  possible. 

We  want  to  evaluate  how  well  the  approximation  of  the  availability  of  the  component  is.  A 
simulation  in  ProModel  simulation  software  is  built  to  compare  the  approximation  to.  The 
modeled  system  seems  simple  enough  to  simulate.  In  thinking  about  the  exact  timing  of 
every  event  in  the  simulation  however,  I  realized  that  discrete  event  simulation  software  is 
not  always  specifically  designed  to  simulate  analytical  models.  Debugging  the  simulation 
with  the  trace-function  is  necessary  to  ensure  that  the  simulation  is  operating  exactly  as  stated 
in  the  model  formulation.  To  obtain  a  5%  confidence  interval  on  our  simulation  data,  we 
calculate  that  around  40.000  replications  are  needed.  The  times  of  failure  of  the  component 
are  written  to  a  text-file,  that  is  imported  in  excel.  [Appendix] 

In  looking  at  the  resulting  graphs  for  the  expected  time  of  failure  of  the  component  we  can 
conclude  that  the  cumulative  derivation  method  is  fairly  accurate.  We  also  see  that  it’s 
accuracy  increases  as  the  truncation  level  increases,  which  seems  to  be  a  logical  observation. 
The  mathematical  tractability  as  the  system  size  increases  is  this  methods  greatest  strength. 
This  semester  is  coming  to  a  rapid  close,  and  I  regret  not  being  able  to  explore  the  advantages 
this  method  has  when  we  increase  the  size  of  the  system.  The  goal  to  expand  the  model  to 
phase-type  distributions  has  also  not  been  realized.  I  realize  that  knowledge  is  the  most 
abundant  resource  in  nature,  and  that  there  is  still  a  lot  of  work  to  be  done.  I  would  like  to 
thank  everybody  who  has  made  it  possible  for  me  to  be  exposed  to  the  practice  of  doing 
research  and  the  field  of  stochastic  modeling  in  particular.  My  interest  is  to  continue  as  a 
graduate  student  in  the  field  of  Industrial  Engineering,  specifically  in  the  areas  of  operations 
research  and  simulation  modeling.  Special  thanks  to  Dr.T.I.  Matis  who  challenged,  inspired, 
and  supported  me  throughout  this  entire  experience. 
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Library  Research 


Cassady,  C.R.,  Iyoob,  I.M.,  “Analysis  of  Equipment  Availability  under  Varying  Corrective 
Maintenance  Models”.  Work  Review  (2002) 

In  this  paper,  the  authors  present  availability  measures  of  repairable  equipment. 
Different  corrective  maintenance  models  are  compared  and  the  effects  of  varying 
parameters  of  the  models  are  examined.  The  effect  of  changing  corrective  maintenance 
models  on  equipment  availability  is  analyzed.  The  point  availability  is  the  major  measure 
of  performance,  which  is  illustrated  in  the  paper  by  various  graphics.  Results  are  obtained 
using  a  simulation  model  created  in  Visual  Basic,  and  a  95%  confidence  interval  is 
obtained  by  multiple  repetitions.  The  models  of  the  impact  of  repair  include  Perfect  Repair, 
Minimal  Repair,  Kijima  Types  Repair,  and  Quasi-Renewal  Repair.  The  simulation  model  is 
developed  to  obtain  functional  approximations  to  the  availability  function  for  each  of  these 
models.  This  research  does  not  provide  any  analytical  solution  methods,  but  it  can  be  used 
as  a  tool  when  measuring  the  relative  error  of  alternative  analytical  methods. 


Dieulle,  Laurence,  “Reliability  of  several  component  sets  with  inspections  at  random  times” 
European  Journal  of  Operational  Research  139  (2002)  96-1 14 

This  paper  considers  a  random  process  representing  a  system  of  components  with 
constant  failure  rates  and  subjected  to  inspections  at  times  defining  a  renewal  process.  An 
analytical  method  for  calculating  the  reliability  function,  its  Laplace  transform  and  the 
mean  time  to  failure  is  given.  These  formulas  are  only  computable  if  the  Laplace  transform 
of  the  inter-arrival  law  of  the  renewal  process  is  explicit.  The  asymptotic  behavior  of  the 
reliability  and  the  failure  rate  of  the  system  are  studied.  The  reduced  ability  to  model  a 
variety  of  distributions,  and  the  intractability  when  modeling  large-scale  systems  however 
makes  this  study  less  applicable  in  reality. 


Kijima,  M.  “Some  results  for  repairable  systems  with  general  repair”  Journal  of  Applied 
Probability  26  (1989)  89-102 

In  this  paper  Kijima  develops  general  repair  models  for  a  repairable  system  by  using 
the  idea  of  the  virtual  age  process  of  the  system.  Two  models  are  constructed  depending  on 
how  the  repair  affects  the  virtual  age  process.  These  models  are  then  used  to  obtain  an 
upper  bound  for  the  expected  value  of  the  survival  function  when  a  general  repair  is  used. 
The  introduction  of  the  Virtual  Aging  Process  is  interesting  especially  considering  the  large 
number  of  references  that  have  been  made  to  this  paper.  Some  knowledge  of  the  survival 
function  makes  modeling  maintenance  policies  easier,  and  the  idea  of  virtual  aging  can 
greatly  contribute  to  this  endeavor. 


Lam,  Y.,  Tony,  H,  K.,  “A  general  model  for  consecutive-k-out-of-n:  F  repairable  system  with 
exponential  distribution  and  (k-l)-step  Markov  dependence.”  European  Journal  of 
Operational  Research  129  (2001)  663-682 

In  this  paper,  a  general  model  for  a  repairable  system  in  which  the  lifetime  of  a 
component  depends  on  the  number  of  consecutive  failed  components  that  precede  the 
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component.  The  failure  and  repair  time  are  both  exponentially  distributed  in  this  model.  A 
transition  density  matrix  is  determined,  and  a  general  Markovian  model  is  used  to  obtain 
some  measures  of  performance  such  as  the  availability,  the  rate  of  failures  and  the 
reliability.  The  assumption  is  that  a  failed  component  after  repair  will  be  “as  good  as  new”. 
This  publication  has  some  general  ideas  that  crossover  to  the  other  references  presented 
here,  and  to  my  own  field  of  research.  The  proposed  methods  however  lack  the  ability  to 
solve  systems  that  show  failure  and  repair  times  that  are  not  exponentially  distributed.  It 
exemplifies  that  there  is  a  serious  need  for  analytical  methods  that  can  provide  measures  of 
performance  for  larger  scale  systems  with  non-exponentially  distributed  random  variables. 


Matis,  J.H.,  Kiffe,  T.R.,  “Stochastic  Population  Models,  a  compartmental  perspective” 

A  classic  book  of  Jacquez  on  compartmental  analysis  inspired  this  publication  in 
which  stochastic  compartmental  analysis  is  reviewed  and  generating  functions  are  used  to 
obtain  many  of  the  results.  The  theoretical  development  of  the  methods  used  is  found  in 
Chapters  three  and  nine,  for  all  practical  purposes  I  focused  on  chapter  three.  This  chapter 
explains  the  use  of  moments  and  cumulants  when  describing  single  population  stochastic 
models.  A  standard  approach  for  solving  probability  functions  known  as  the  Kolmogorov 
differential  equations  is  illustrated.  Cumulant  functions  are  shown  to  be  very  useful  for 
finding  a  distributions  cumulants  by  obtaining  and  solving  partial  differential  equations  for 
associated  generating  functions.  An  immigration  death  model  is  used  to  illustrate  the 
theoretical  development  of  the  generating  function  approach.  The  use  of  cumulant 
functions  presented  in  this  analysis  is  applied  in  the  arena  of  reliability  by  Matis  and 
Feldman,  as  described  above.  I  have  extended  these  ideas  to  find  practical  application  in 
reliability  systems  with  the  support  of  Dr.  T.I. Matis. 


Matis,  T.I.,  Feldman,  R.M.,  “Using  Cumulant  Functions  in  Queueing  Theory”  Queueing 
Systems  40,  (2002)  341-353 

This  publication  demonstrates  a  new  procedure  for  obtaining  measures  of 
performance  of  state-dependent  queueing  networks.  This  procedure  uses  cumulant 
generating  functions  and  relates  these  to  the  intensity  functions  in  the  network.  The  service 
rate  is  expressed  as  polynomial  function  of  the  state  of  the  system,  from  which  a  partial 
differential  equation  of  the  cumulative  generating  function  is  obtained.  This  partial 
differential  equation  then  yields  a  set  of  ordinary  differential  equations,  which  are  then 
solved  to  obtain  the  first  and  second  moments  of  the  system  with  Markovian  arrival  and 
service  rates.  The  first  and  second  moments  of  the  random  variables  in  the  system  provide 
transient  information  when  describing  the  system.  As  the  network  increases  in  size  this 
solution  method  remains  tractable,  in  contrast  to  the  method  proposed  by  Rehmert  in  the 
dissertation  reviewed  above.  Cumulant  functions  have  previously  been  used  by  Matis  and 
Kiffe  to  obtain  measures  of  performance  of  ecological  models,  such  as  the  spread  of 
honeybees  or  muskrats.  This  cumulant  derivation  procedure  is  also  used  in  the  research  of 
reliability  systems  subject  to  non-fatal  shocks. 
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Pham,  H.  Wang,  H.  “Imperfect  maintenance”  European  Journal  of  Operational  Research  94 
(1996)  425-438 

The  maintenance  of  deteriorating  systems  is  often  imperfect,  as  studied  by  using 
simulation  in  the  first  reference  by  Cassady  and  Iyoob.  Imperfect  maintenance  studies 
using  mathematical  models  for  estimating  availability  functions  and  reliability  have 
undergone  some  breakthroughs  that  are  discussed  and  summarized  in  this  paper.  Treatment 
methods  for  imperfect  maintenance,  such  as  the  (p,  q)  and  (p(t),  q(t))  rules,  the 
improvement  factor  method,  virtual  age  method,  shock  model  method,  (□  □)  rule,  are 
examined  and  explained.  Preventative  maintenance  policies,  such  as  age-dependent, 
periodic  PM,  failures  limit  policy,  sequential  PM  policy,  repair  limit  policy,  and  multi- 
component  systems  can  indicate  what  model  needs  to  be  selected.  Although  the  focus  of 
this  paper  tends  towards  modeling  maintenance  policies  to  reduce  cost,  the  underlying 
principals  can  be  very  useful  when  thinking  about  reliability  systems  in  a  more  general 
sense. 


Rehmert,  I.  J.,  “Time-Dependent  Availability  Analysis  for  the  Quasi-Renewal  Process.”  Diss. 

Virginia  Polytechnic  Institute  and  University  (2000) 

This  research  is  based  upon  the  quasi-renewal  process  proposed  by  Wang  and  Pham 
which  is  an  alternative  to  the  widely  studied  imperfect  repair  model,  the  (p,  q)  model, 
proposed  by  Brown  and  Proschan.  A  quasi  renewal  process  can  realistically  describe  the 
behavior  of  repairable  equipment.  The  framework  provided  in  this  dissertation  allows  for 
the  description  of  the  time-dependent  behavior  of  this  non-homogeneous  process.  Two 
equivalent  expressions  for  the  point  availability  of  a  system  with  operation  intervals  and 
repair  intervals  that  deteriorate  according  to  a  quasi-renewal  process  are  constructed.  These 
expressions  are  used  to  provide  upper  and  lower  bounds  on  the  approximated  point 
availability.  Laplace  transforms  are  used  to  solve  the  resulting  expressions.  The  quasi 
renewal  function  and  the  point  availability  function  are  found  for  exponential,  normal,  and 
gamma  operating  and  repair  intervals.  It  is  necessary  to  truncate  the  expressions  to  invert  to 
the  time  domain  to  obtain  numerical  results.  The  usefulness  of  this  approach  seems  limited 
because  it  does  not  find  exact  expressions  for  all  distributions,  and  the  calculations  seem  to 
be  intractable  as  the  size  of  the  system  increases. 


Scarsini,  M.  Shaked,  M.  “On  the  value  of  an  item  subject  to  general  repair  or  maintenance” 
European  Journal  of  Operational  Research,  Vol.122  (2000)  625-637 

This  paper  introduces  and  studies  finding  a  practical  expression  of  the  monetary  value 
of  an  item.  The  model  that  is  constructed  takes  the  repair  and  the  maintenance  procedures 
that  are  applied  to  it  during  its  lifetime.  The  number  of  repairs  and  the  degree  of  the  repair 
are  taken  into  account  when  doing  the  calculations.  The  degree  of  the  repair  follows  the 
ideas  of  virtual  aging  modeling  as  was  presented  earlier  by  Kijima  in  his  virtual  aging 
models.  The  uncertainty  is  introduced  into  the  model  by  the  distributions  of  the  inter  repair 
or  inter  maintenance  periods.  This  paper  shows  how  the  virtual  aging  process  can  be  used 
to  model  cost  under  repair  and  maintenance  procedures. 
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Sheu,  S.,  Griffith,  W.S.,  “Multivariate  Age-Dependent  Imperfect  Repair”  Naval  Research 
Logistics,  Vol.  38  (1991)  839-850 

This  article  considers  models  of  systems  whose  components  have  dependent  life 
lengths  with  specific  multivariate  distributions.  Components  are  repaired  according  to  a 
corrective  maintenance  scheme,  meaning  that  components  are  repaired  upon  failure.  Only 
two  types  of  repair  are  considered  in  this  paper:  perfect  repair,  and  imperfect  repair.  The 
study  focuses  on  a  model  in  which  the  nature  of  the  repair  is  age  dependent.  The  model 
uses  the  (p(t),  q(t))  rule  which  was  earlier  described  in  this  reference  list  in  the  publication 
by  Wang.  An  expression  for  the  cumulative  hazard  function  is  derived,  which  can  be  useful 
when  describing  reliability  of  systems.  The  paper  however  lacks  any  numerical  examples 
that  demonstrates  finding  the  hazard  function. 


Wang,  H.,  “A  survey  of  maintenance  policies  of  deteriorating  systems”  European  Journal  of 
Operational  Research  139  (2002)  469-489 

This  survey  summarizes,  classifies,  and  compares  various  existing  maintenance 
policies  for  both  single-node  and  multi-node  systems.  All  these  models  fall  into  categories 
such  as:  Age  replacement  policy,  block  replacement  policy,  periodic  preventive 
maintenance  policy,  failure  limit  policy,  sequential  preventive  maintenance  policy,  repair 
cost  limit  policy,  repair  time  limit  policy,  repair  number  counting  policy,  reference  time 
policy,  mixed  age  policy,  group  maintenance  policy,  etc.  The  characteristics,  advantages, 
and  drawbacks  for  each  kind  of  policy  are  addressed.  Maintenance  and  replacement 
problems  have  been  studied  for  the  past  several  decades,  and  this  invited  review  gives  a 
clear  overview  of  the  models  used.  Wang’s  research  has  introduced  the  concept  of  “Virtual 
Age”  when  modeling  Quasi-renewal  or  imperfect  repair. 
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Based  on  9990  repetitions.  Failure  rate  is  related  to  both  current  and  cumulant  count 


CD 

O 

o 

\0 

O'' 

00 

NP 

O' 

CO 

CM 

o 

in 

in 

o 

CO 

CM 

CM 

T— 

O 

fs. 

vP 

o' 

sO 

O'* 

CD 

CM 

CO 

in 

o 

CO 

CD 

CO 

CO 

CM 

CM 

«0  O  g  t 

rqogo 

II  'U  0, 1  ® 
a.  a.  S  oq. 


>1  & 

W  UJ 


000 


000 


00.  00 


00.  00 

Scenario  :  Normal  Run 

Replication  :  999  of  999 
Simulation  Time  :  0.002319444444  hr 


LOCATIONS 


Location  Scheduled 
Name  Hours 

Environment  0.002319444444 
Machine  0.002319444444 
Current  0.002319444444 
Choose  0.002319444444 
Write  up  0.002319444444 


Average 
Total  Seconds 
Capacity  Entries  Per  Entry 


999999  2  0.000000 
999999  3  1.226667 
999999  2  1.650000 
999999  1  0.000000 
999999  1  0.000000 


Average 

Maximum 

Current 

Contents 

Contents 

contents 

0 

1 

0 

0.440719 

3 

3 

0.39521 

2 

1 

0 

1 

1 

0 

1 
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LOCATION  STATES  BY  PERCENTAGE  (Multiple  Capacity) 


Location 

Name 

Scheduled 

Hours 

% 

Empty 

/o 

Partially 

occupied 

% 

Full 

% 

Down 

Envi ronment 

0.002319444444 

100.00 

0.00 

0.00 

0  00 

Machine 

0.002319444444 

64.91 

35.09 

0.00 

0.00 

Current 

0.002319444444 

64.91 

35.09 

0.00 

0.00 

Choose 

0.002319444444 

100.00 

0.00 

0.00 

0.00 

write  up 

0.002319444444 

100.00 

0.00 

0.00 

0.00 

RESOURCES 


Resource 

Name 

Units 

scheduled 

Hours 

Number 
of  Times 
Used 

Average 

seconds 

Per 

Usage 

%  util 

Repai r . 1 

i 

0.002319444444 

1 

2.930000 

35.09 

Repai r . 2 

i 

0.002319444444 

1 

0.370000 

4.43 

Repai r . 3 

i 

0.002319444444 

0 

0.000000 

0.00 

Repair. 4 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r . 5 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r. 6 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r .7 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r . 8 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r . 9 

i 

0.002319444444 

0 

0.000000 

0.00 

Repai r . 10 

i 

0.002319444444 

0 

0.000000 

0.00 

to 
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«  momcum .  m 

UptoOrder  [nVar_,  ord._J  :  = 

Rest  [Select  [Distribute  [Table  [Range  [0,  ord]  ,  {nVar}]  ,  List]  ,  (Plus  @@  #  ^  ord)  &]  ]  ; 
orderListl =  UptoOrder [3,  3]  ? 

relations  mod  = 

MomCumConvert  [#,  ForMomentQ  -4  "Y"  ,  CenteredQ  ->  "N"  ,  HomentSymbol  M,  CumulantSymbol  -4  K]  &  /@ 
orderListl; 

relations  :  =  relationsmod  /  .  {K[i _ ]  -4  K[i]  [t]  } 

momcumrule  =  relations  /  .  {Equal  -4  Rule)  ; 
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:=h/,  {©i  -►  Evaluate  [i*ax]  ,  ©2  -4  Evaluate^  *a2]  ,  ©3  -*  Evaluate  [k  *  ©3  ]}  ; 
duml  [i_ ,  j_ ,  k_J  :  =  Coefficient  [Evaluate  [ 

eql[Sign[i],  Signfj]  ,  Sign[k]]],  0102  03]  / .  momcumrule  ; 
eq2  [i_,  j_,  3c_]  :  =  g  /  .  {0i  ->  Evaluate[i  *  ©i]  ,  ©2  Evaluate  [j  *  ©2  ]  ,  ©3  -♦  Evaluate  [k  *  ©3]  }  ; 
dum2[i_,  j_,  k_]  :=  Coefficient  [Evaluate  [ 

eq2 [Sign [i] ,  Sign[j] ,  Sign[k]]],  e{  e\  ©3]  /  .  momcumrule; 

bilbo  =  Map  Thread  [duml.  Transpose  [orderListl]  ]  ; 

bilbo2  =  MapThread  [dum2 ,  Transpose  [orderListl]  ]  ; 

Shi  take  =  Table  [bilbo2  [  [i]  ]  ==  bilbo  [  [i]  ]  ,  {i,  1,  19}]  ; 

Bonzai  =  Tablet 

K [Part  [orderListl,  i,  1]  ,  Part  [orderListl,  i,  2]  ,  Part  [orderListl,  i,  3]]  '  [t]  ,  {i,  1,  19}]  ; 
neweqns  =  Solve  [Shi take,  Bonzai]  ; 
neweqmod  =  neweqns  /  .  {Rule  ->  Equal}  ; 

Samuri  =  Table  [K  [Part  [orderListl,  i,  1]  ,  Part  [orderListl ,  i,  2]  ,  Part  [orderListl,  i,  3]]  [0]  =  0, 
{±,1,19}]; 

Monkey  =  ReplacePart [Samuri,  K[0,  0,  1]  [0]  ==  1,  1]  ; 

Joy  =  Join [First [neweqmod] ,  Monkey] ; 

Pokemon = 

Table  [K  [Part  [orderListl,  i,  1],  Part  [orderListl,  i,  2],  Part  [orderListl ,  i,  3]],  {i,  1,  19}]; 

rs  =  NDSolve[  Joy ,  Pokemon,  {t,  0,  50}  ,  MaxSteps  -4  10000] 

{{K[0,  0,  1]  -4 InterpolatingFunction[ { {0 . ,  50.}},  <>], 

K[  0 ,  0,  2]  *4  Int erpolatingFunction [ { { 0 . ,  50.}},  <>], 

K [ 0 ,  0,  3]  -4 InterpolatingFunctionf { {0. ,  50.}},  <>], 

K[0,  1,  0]  Int erpol at ingFunct i on [ { { 0 .  ,  50.}},  <>], 

K [ 0 ,  1,  1]  -4 InterpolatingFunctionf {{ 0 . ,  50.}},  <>], 

K[0,  1,  2]  <4  Int  erpol  at  ingFunct  ion  [{ {0 . ,  50.}},  <>], 

KfO,  2,  0]  -4 InterpolatingFunctionf {{ 0 . ,  50.}},  <>] , 

K[0,  2,  1]  *4  Int  erpol  at  ingFuncti  on  [  {  { 0  .  ,  50.}},  <>]  , 

K[0,  3,  0]  ->  InterpolatingFunctionf  { {0. ,  50.}},  <>], 

K  [  1 ,  0,  0]  -4  InterpolatingFunctionf  {{ 0 . ,  50.}},  <>], 

K[  1 ,  0,  1]  Int erpol at ingFuncti on [ { { 0 . ,  50.}},  <>], 

K  [  1 ,  0,  2]  -4  InterpolatingFunctionf  {{0.  ,  50.}},  <>], 

Kfl,  1,  0]  -> InterpolatingFunctionf { (0 . ,  50.}},  <>], 

Kfl,  1,  1]  -4  InterpolatingFunctionf { {0. ,  50.}},  <>], 

Kfl,  2,  0]  -4  InterpolatingFunctionf {{ 0 .  ,  50.}},  <>], 

K [ 2 ,  0,  0]  -4  InterpolatingFunctionf { {0 . ,  50.}},  <>]  , 

K[2,  0,  1]  -4  InterpolatingFunctionf { {0. ,  50.}},  <>], 

K  f  2 ,  1,  0]  -4  InterpolatingFunctionf  {  {0  . ,  50.}},  <>], 

K[3,  0,  0]  -4  InterpolatingFunctionf { {0 . ,  50.}},  <>]}} 
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K [0 ,  0,  1] [2]  /.  rs 
K[0,  0,  1]  [4]  /.  rs 
K[0,  0,  1]  [6]  /.  rs 
K [0 ,  0,  1] [8]  /.  rs 
K[0,  0,  1J  [10]  /.  rs 
K[0 ,  0,  1]  [12]  /.  rs 
K [0 ,  0,  1]  [14]  /.  rs 
K [0 ,  0,  1]  [16]  /.  rs 
K [0 ,  0,  1]  [18]  /.  rs 
K [0 ,  0,  1]  [20]  /.  rs 
K [0 ,  0,  1] [22]  /.  rs 
K [0 ,  0,  1]  [24]  /.  rs 
K[0,  0,  1]  [26]  /.  rs 

Plot  [Evaluate  [{K[l,  0,  0]  [t]  }  /.  rs]  ,  {t,  0,  50}  , 
PlotRange  -♦  {0 ,  10}  ,  AxesLabel  -»  {  "  t"  ,  "E[Xi  (t)  ]  "  }  ] 
Plot  [Evaluate  [{K[0,  1,  0]  [t]  }  /.  rs]  ,  {t,  0,  50}, 
PlotRange  -►  {0 ,  5}  ,  AxesLabel  -»  {"  t"  ,  "E  [X2  (t)  ]  "  }  ] 
Plot  [Evaluate  [{K[0,  0,  1]  [t]  }  /.  rs]  ,  (t,  0,  50}  , 
PlotRange  ->  {0 ,  1}  ,  AxesLabel  -♦  { "  t"  ,  "E  [X3  (t)  ]  "  }  ] 
Plot  [Evaluate  [  {K  [2 ,  0,  0][t]}  /.  rs]  ,  {t,  0,  50}, 
PlotRange  ->  {0 ,  10}  ,  AxesLabel  -»  { "  t"  ,  "  Var  [Xi  (t)  ]  "  }  ] 
Plot  [Evaluate  [{K[0,  2,  0]  [  t]  }  /.  rs]  ,  {t,  0,  50}  , 
PlotRange  -»  {0 ,  5}  ,  AxesLabel  -+  { "  t"  ,  "Var  [X2  (t)  ]  "  }  ] 

{0.966056} 

{0.858759} 

{0.667377} 

{0.430791} 

{0.219442} 

{0.0838102} 

{0.0228014} 

{0.00419834} 

{0.000497058} 

{0.0000359521} 

{1.50969x10-®} 

{3.50847x  10-®} 

{4.46547  x  lO-10 } 
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Abstract 


The  reliability  modeling  of  numerous  physical  systems  is  critical  in  the  prevention  of  system  failure.  In 
many  instances,  the  failure  rate  of  system  components  is  a  function  of  non-fatal  shocks  or  stresses  to  the 
system  that  occur  at  discrete  points  in  time.  These  shocks  are  assumed  to  be  identical  and  reparable,  and 
they  impact  the  failure  rate  of  the  system  in  a  non-linear  fashion  via  the  cumulative  and  current  number  of 
shocks.  In  this  paper,  we  demonstrate  the  application  of  the  cumulant  derivation  procedure  to  this 
reliability  system  in  a  Markovian  environment  This  approach  utilizes  a  truncated  cumulant  generating 
function  to  generate  a  set  of  ordinary  differential  equations  whose  numerical  solution  approximates  the 
reliability  function.  These  approximations  are  obtained  under  various  truncation  levels  whereby  this 
approach  is  shown  to  be  tractable  for  large  systems. 

1.  Introduction 

In  this  paper,  we  consider  a  single  component  operating  in  an  environment  of  non-fatal  reparable 
shocks.  These  shocks  may  represent  a  wide  variety  of  events,  including  the  failure  of  other  system 
components,  instantaneous  system  stresses,  or  the  states  of  a  compartmentalized  wear  process.  The  shocks 
are  assumed  to  be  homogenous  and  their  arrival  is  governed  by  a  stationary  Poisson  process.  They  are 
repaired  according  to  an  exponentially  distributed  infinite  server  repair  process  that  begins  immediately 
upon  shock  arrival  whose  rate  is  a  function  of  the  total  number  of  cumulative  shocks.  The  failure  of  the 
component  is  Poisson  distributed  with  the  rate  being  a  function  of  both  the  current  number  of  unrepaired 
shocks  and  the  total  number  of  shocks  that  have  ever  been  received  by  the  system.  In  other  words,  each 
shock  has  both  an  immediate  reparable  effect  and  a  permanent  weakening  effect  on  the  system.  Let  X]  (t) 
and  X2(t)  be  integer-valued  random  variables  taking  values  in  [0,1,.  .,oo]  that  denote  the  cumulative  and 
current  number  of  shocks  at  time  t  respectively.  Let  X3(t)  be  an  integer-valued  random  variable  taking 
values  [0,1]  that  represents  the  state  of  the  component  at  time  t,  where  X3(t)=l  denotes  a  functioning 
component  and  X3(t)  =0  denotes  a  failed  component.  All  possible  unit  changes  that  may  occur  in  the  state 
of  the  system  in  a  small  interval  of  time  are  contained  in  the  set  B  and  the  corresponding  state-dependent 

intensity  (rate)  functions  will  be  denoted  as  f(bi,b2,b3)(X,,X2,X3)  for  (b,,b2,b3)eB  .  This  system  is  graphically 
depicted  in  Figure  1. 

We  assume  that  initially  the  component  is  operational,  X3(0)=1,  and  no  shocks  have  been  received 
Xi(0)=X2(0)=0.  The  intensity  function  corresponding  to  component  failure,  f(0,o,-i)(X,,X2,X3),  is  specified 
as  a  nonlinear  function  of  the  cumulative  and  current  number  of  shocks,  X,  (t)  and  X2(t),  and  represents  the 
damage  process  previously  described.  Our  primary  interest  is  in  approximating  the  reliability  function  of 
the  component,  i.e.  the  expected  value  of  X3(t),  for  all  t  ^  0  using  truncated  cumulant  generating  functions. 
These  approximations  are  compared  to  simulated  values  for  several  systems  under  various  intensity 
function  specifications  and  truncation  levels. 


(1) 


Figure  1:  Graphical  Representation  of  Reliability  System 


2.  Cumulant  Derivation  Procedures 

Let  X(t)=(X,  (t),X2(t)„X3(t),)  be  a  random  vector  of  the  system  state  at  time  t.  It  follows  that  X(t) 
forms  a  Markov  process  with  an  absoibing  state  that  denotes  component  failure.  While  such  a  process  may 
e  solved  exactly  using  convention  means,  i.e.  Kolomogorov  Equations,  such  an  approach  is  intractable  for 
-ff?.  ®ystemf  or  111086  with  311  infinite  state  space.  As  an  alternative,  the  cumulants  of  the  state-distribution 
of  X(t)  may  be  approximated  using  a  truncated  generating  functions.  These  cumulant  measures  correspond 
directly  to  the  common  measure  of  mean,  variance,  covariance,  skewness,  etc.  of  the  state  distribution. 
Previous  investigations  into  the  nature  of  cumulant  functions  by  Kendall[l][2]  and  Smith[3]  reveal  several 
interesting  properties.  In  particular,  the  cumulants  of  a  multivariate  normal  distribution  greater  than  the 
second  order  are  null,  and  the  marginal  cumulants  of  a  Poisson  distribution  are  equal.  This  and  other 
properties  of  cumulants  are  exploited  by  this  approximation  procedure.  This  section  contains  only  a  brief 

°f  ^  cumulant  derivation  procedure  based  on  the  full  development  found  in  Matis  and  Feldman 
[4][5],  Matis  [6],  and  Matis  and  Kiffe  [7], 

Let  M(6X  ,02,03 ,t)  be  the  multivariate  moment  generation  function  of  X(t)  defined  in  the  usual 
manner,  and  let  K(0u02,03,t)  be  the  multivariate  cumulant  generating  function  defined  as 


K(0v02,0vt)=  £ 

#1  >• 


kav...,anml-<n 


where  A,r+  denoted  the  set  of  non-negative  integers.  The  joint  cumulants  a  (?)  of  the  system  are 

defined  as  functions  of  the  individual  moments  through  the  relationship  between  the  generating  functions, 
1.6. 


The  moment  generating  function  of  X(t)  related  to  the  polynomial  intensity  functions  of  the 
system  through  a  partial  differential  equatioa  This  relationship  was  investigated  by  Bartlett[8]  and 
Bailey[9]  and  thereby  dubbed  the  “Random  Variable  Technique”  While  the  specification  of  the  partial 
differential  equation  of  the  moment  generating  is  almost  immediate,  finding  a  solution  is  usually 
computationally  intractable.  The  cumulant  derivation  procedure  involves  substituting  a  truncated  cumulant 
generating  function  for  the  moment  generating  function  into  the  partial  differential  equation  according  to 
Eq.  (1).  A  set  of  approximating  ordinary  differential  equations  is  obtained  upon  expanding  the  partial 
derivatives,  substituting  Taylor  series  expansions  for  the  exponential  terms,  and  equating  the  coefficients  of 


unique  combinations  of  the  02  s.  The  accuracy  of  these  approximations  is  dependent  upon  the 

specification  and  order  of  the  polynomial  intensity  functions  and  the  level  of  cumulant  truncation  Previous 
investigations,  Matis[6],  have  shown  that  truncation  at  the  3rd  order  is  generally  sufficient  for  “good” 
approximations  of  the  first  order  cumulants  (mean),  while  that  at  the  4S  order  is  sufficient  for  the  marginal 
second  order  cumulants  (variance).  Further  developments  of  the  cumulant  derivation  procedure  will  be 
descnbed  in  a  reliability  context. 

3.  Cumulant-Based  Analysis  of  the  Non-Fatal  Shock  Process 

fn  this  section,  we  demonstrate  the  application  of  cumulant-based  procedures  to  the  non-fatal 
shock  process  descnbed  in  the  introduction  of  this  paper,  see  Figure  1.  This  will  be  shown  for  one  instance 
of  the  problem  under  a  unique  set  of  intensity  functions.  Let  the  intensity  function  of  the  system  be  defined 


f(i,i,o)(Xi,X2,X3)  =  X 
f:o,-i.o)(X1,X2,X3)  =  p2X,2(t)X2(t) 
f(o,o,.)(X,,X2,X3)  =  p3(X3(t)+X]  (t)X22(t)  X3(t)) 

X=,10’  ^r  05’ afd  t13./025'  111 0ther  words’  die  rate  of  repair  is  dependent  upon  the  cumulative 
number  of  shocks  and  the  failure  rate  of  the  component  is  dependent  on  both  the  current  and  cumulative 

S  ’  <  r  m  manner-  A  partial  differential  equation  of  the  moment  generating 

function  of  X(t)  is  found  using  the  "Random  Variable  Technique”  as  &  & 

mevevevo  =  +e2-\)  ^  _1}  dzM{eveveyt) 

&  d6\  ^  ddfddj 

+  e(-03  -l)f dM{6v0veyt)  dAM{ev62,0yt)\  (2) 

M3  K  d03  d0l8022d03  y 

An  m*  order  truncated  cumulant  generating  function  and  a  Taylor  series  expansion  of 

the  exponential  terms  is  substituted  into  Eq.(2)  yielding  the  expression 


deK^M!'>  =  ( A  (-  e2 )  ) 

&  /!  J  d6\  I  ^d02 


f A  (-  0 3  )  T  deKMM't'>  aV^.*2,V) 

2j—t-  — — +— — ^ — 


V>=1  '•  X  °°3  ddl802L803  J 

Expanding  the  partial  derivatives,  converting  moments  to  cumulants  (Smith[3J),  and  equating  the 
coefficients  of  like  polynomial  terms  on  the  right  and  left  hand  side  of  Eq.  (3)  yields  a  closed  set  of 

r  1  3  \ 

3[  IT  ^  +  m)  I  - 1  ordmaiy  differential  equations.  The  size  of  the  generated  sets  of  ordinary  differential 


equation  does  not  permit  their  demonstration  in  this  paper,  yet  the  number  of  such  equation  is  9  under  a 
truncation  level  of  m-2, 1 9  under  w=3,  and  34  under  m= 4.  Approximations  to  the  low  order  cumulants  are 
found  upon  numerically  solving  these  sets  of  ordinary  differential  equations. 

4.  Numerical  Results 

The  sets  of  ordinary  differential  equations  generated  from  Eq.  (3)  were  numerically  solved  using  the 
mathematical  software  Mathematica®  under  the  truncation  levels  m=2,  3,  4.  These  were  then  compared  to 
simulated  point  estimates  based  on  10,000  replications  using  the  software  ProModel®.  The  initial 


conditions  for  the  process  have  all  cumulants  set  equal  to  zero  except  for  ko,o,i(0),  which  corresponds 
directly  to  the  reliability  of  the  component  at  time  0,  is  set  to  one.  As  previously  noted,  truncating  the 
cumulants  at  m=2  implies  that  the  state-distribution  is  normally  distributed.  Increasing  the  level  of 
truncation  to  m=  3  brings  in  skewness  and  m= 4  brings  in  kurtosis,  in  addition  to  die  effects  of  the  higher- 
order  cross  cumulants.  The  cumulant  approximations  for  the  reliability  of  the  component,  i.e. 
R(t)=E[X3(t)],  is  given  in  Figure  2. 


Figure  2:  Graph  of  R(t)  =  E[X3(t)]  for  Varying  Truncation  Levels 


The  approximation  of  R(t)=E[X3(t)]  is  relatively  tight  between  all  values  of  m  and  close  to  the  simulated 
value.  This  result  is  promising  as  the  reliability  function  for  systems  subject  to  similar  non-fatal  shock 
processes  may  be  well  approximated  using  small  values  of  m.  As  such,  this  result  provides  evidence  that 
approach  may  be  extended  to  similar  large-scale  networks  in  a  computationally  efficient  manner.  Though 
not  of  primary  interest  in  this  paper,  the  differences  in  approximations  of  the  variance  of  the  current 
number  of  shocks,  Var[X2(t)],  between  truncation  levels  is  noteworthy  and  is  given  in  Figure  3. 


Var[X2(t)J 


Figure  3:  Graph  of  Var[X2(t)]  for  Varying  Truncation  Levels 


The  specified  model  was  then  evaluated  under  the  parameter  specification  X=.  10,  p2  =  05,  and  p3 
=.005  yielding  the  graph  in  Figure  4.  Comparing  the  similarity  of  the  cumulant  approximations  for  various 
truncation  levels  in  both  Figures  2  and  4  provides  evidence  that  the  precision  of  the  reliability  function 
approximations  do  not  significantly  vary  with  the  rate  of  system  failure.  Comparing  these  approximations 
to  simulations,  however,  provides  evidence  that  the  accuracy  of  the  approximations  increases  as  the  rate  of 
system  failure  decreases. 


E[Xj(t)] 
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Figure  4:  Graph  of  R(t)  =  E[X3(t)]  for  Vaiying  Truncation  Levels 

The  polynomial  intensity  function  of  the  model  were  then  redefined  as 

fci,i,o)(Xi,X2>X3)  =  X 

%>,-!, 0)(X,,X2,X3)  =  p2X,3(t)X2(t) 

Wn(X,,X2,X3)  =  g3(X3(t)+X,2(t)X23(t)  X3(t)) 


-  m  =  2 
m  =  3 
m  =  4 
■  sim 


increasing  the  dependency  of  the  failure  and  repair  rates  on  the  state  of  the  network.  The  graph  of  the 
reliability  function  is  given  in  Figure  5.  This  deviation  of  the  cumulant  approximation  under  m  2  provides 
evidence  that  2nd  order  truncation  is  not  sufficient  for  approximating  the  cumulants  of  systems  that  are 
strongly  state-dependent.  This  observed  result  is  consistent  with  the  previously  stated  properties  of 
cumulants,  i.e.  truncation  at  m=2  assumes  a  multivariate  normal  state  distribution.  This  normal  assumption 
clearly  does  not  hold  for  systems  with  strong  state-dependency  in  which  skewness  is  clearly  present. 
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Figure  5:  Graph  of  R(t)  =  E[X3(t)]  for  Varying  Truncation  Levels 


5.  Conclusions 

In  this  paper,  we  have  shown  the  application  of  cumulant  derivation  procedures  to  a  reliability 
system  subject  to  non-fatal  shocks,  i.e.  a  state-dependent  reliability  system.  The  effect  of  various 
truncation  levels  on  the  accuracy  of  the  reliability  function  approximation  was  demonstrated  for  various 
parameters  and  intensity  functions.  The  similarities  between  these  approximations  provide  insight  into  the 
expandability  of  the  approach  to  large,  complex  systems.  A  copy  of  the  Mathematica®  computational 
routines  used  to  set  up  and  solve  the  approximating  set  of  ordinary  differential  equations  may  be  obtained 
from  the  authors  upon  request. 
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